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Abstract 

The prediction of extremely high wind speeds , at ground 
level on the downstream side of a mountain range, is possible 
by solving the initial value problem for a two-layered nonlinear 
"shallow-water" model of the atmosphere. Three different 
numerical methods are described to find the solutions which may 
involve shocks: (i) the vonNeumaun-Richtmyer artificial viscosity 

method, (ii) a altering scheme, and (iii) a hybrid method. 
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1. I ntroduction 

Occasionally, when air is crossing at moderate speed 
over the top of a mountain range, extremely high wind speeds 
will develop at ground level on the lee (downstream) side 
(51 . Simple one-layered [4] and two-layered [31 
nonlinear "shallow-water" models of the troposphere were 
designed and used to study this phenomenon. For a large 
class of initial values, such asymmetric steady flows were 
found to develop quickly in time, even when the mountain 
is taken to be a symmetric parabolic arc. The numerical 
method used in [3] to produce tho steady-state solutions 
made use of a von Neumann and Richtmyer [61 pseudo-viscosity 
term, which was "switched on" in regions of "compression" 
only. This artificial viscosity method was needed to suppress 
the post-shock oscillations that would have otherwise 
developed in the levels and velocities of the two-layered 
model. Since these flows often contained stationary shocks 
on the lee side of the mountain, such post shock oscillations 
would have seriously masked the true nature of the steady 
states. In [3J, the use of this "compression switch" 
prevented the application of the artificial viscosity term 
near the upstream base of the mountain, for a physically 
realizable subset of initial values- In these cases, spurious 
oscillations developed on the upstream side of the mountain 
and continued downstream, so that the detection of a steady state 
flow over the mountain was prevented. 
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The present work verifies that- indeed for all initial 
values of physical interest, steady states do develop over 
the mountain in the two-layered model. We therefore believe 
that the shallow-water models can be used to forecast such 
ground level strong winds on the lee side of mountains. 

The correct steady state solutions are found 
independently by the tfge of three different numerical 
methods: (i) by eliminating the compression switch, i.e., 
so that the pseudo-viscosity term is always present; 

(ii) by using a filtering scheme; (iii) by using a hybrid 
scheme. The numerical results obtained from these three methods 
have close agreement In section 7 , we describe the differential 
equations of the two-layered model, the parameters that 
characterize the initial values, and the nature of the 
physically relevant solutions. In section 3, we describe 
the three numerical methods and comment on their advantages 
and disadvantages. 
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2 . Differential equations, initial data, and solutions , 


The time dependent, two-layered, long wave, incompressible 
fluid model, for flow in a plane perpendicular to the earth 
and to the mountain ridge, is governed by the partial differ- 
ential equation system [3] , for the vector w * w(x,t) : 
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We have chosen u and 4 to represent the horizontal velocity and 
the depth respectively of the lower’ fluid layer, and the same 
quantities primed (u • r <4 * ) to represent velocity and depth 
for the upper layer (see Figure 1) . Since we will be computing 
flows that have shocks, it is appropriate to ;/ork with 


m=u$, m* = u*4* » 
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representing the momentum per unit length (and density) , in 
the lower and upper layers respectively. The constant quantity, 
r (r £ 1) , is the ratio Pj/P of the densities of the 
upper layer and p of the lower layer, g denotes the 
gravitational constant, while H = H(x) is the formula represent- 
ing the height, above some datum, of the ground. We remark 
that [3] treats also the case of a three-layered model, 
in which the highest layer is passive, by means of a quite 
similar system of equations that incorporates a constant, 
s; namely, the ratio of the densities of the topmost and bottom- 
most layers and that deals solely with the same dependent 
variables u, <{>, u', and <J> a . It would have been simple to 
formulate the corresponding difference schemes for this 
three-layered model, since it is governed again by a system 
of four, first order, partial differential equations with 
different constants in the coefficients. 

As was observed in [3], the system (1) is hyperbolic, 
when its four characteristic speeds, dx/dt = p , are 
real and distinct; namely, when the roots, p = Pj(x,t) for 
j = 1,2, 3, 4, of the following quartic equation are real and 
distinct, < 2 < y 3 < • 

(2) [(u - p) 2 - g<f>Jl(u' - p) 2 - g$ ’ J - rg 2 <f> = 0 . 

We expect that the initial value problem for (1) is well posed, 
when the four roots are real and distinct for the initial data. 
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The family of steady solutions (i.e., time independent) 
of (1) is found in [3] by studying the limit as t • 

r 

of the numerical solutions of initial value problems for 
suitable choices of impulsive initial data. That is, the 
initial velocities u(x,0) and u*(x<0) are set to be constant 
and the tops of the two layers are initially chosen to be 
horizontal. It was found by doing many subsidiary numerical 
calculations with such data, that all of the qualitative 
features of the steady state solutions were exhibited by the 
solutions that arose from a special choice of the parameters 
that represent the initial data. That is, when the mountain 
profile is a simple parabolic arc 


(3) 

where 



M < a , 

|x| > a , 


(4) 



simply assume that the initial velocities are constant and 
equal , say 


(5) u(x,0) = u'(x,0) = u -co ; 

and suppose that the initially constant height of the lower layer, 
<|>(x,0) + H(x), equals the initially constant depth of the 
upper layer, $*(x,0): 


(6) 


♦ (x,0) + H (x) •* <p 1 (x, 0) « j 
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and finally select a value for r. It then is appropriate 

to consider the two dimensionless parameters 

< 

i 

u - H „ 

(7) F 0 = : , and M c = — j£- 

as determining the initial data. 

In [3],, the plane is divided into regions, 

each producing a different kind of limiting solution as t + 
from the corresponding choice of initial data. In particular, 

we reproduce a figure of 13] as Figure 2, and observe 

* 

that r » 0.8 and regions A, B and B' are the ones that 
correspond to physically realizable velocities in the atmos- 
sphere. We note that this model should not be used to make 
calculations for data in region A. For region A, it would be 
appropriate to use the more conventional linearized theory 
in two spatial dimensions, to describe the sinusoidal lee waves 
that occur in nature. On the other hand, it is the development 
of a high lee side velocity in the steady asymmetric solutions 
for regions B and B', that leads us to recommend the use of 
this nonlinear model to describe and predict the strong wind 
phenomena on the downstream side of a mountain range. 

In [3] , a steady solution was approached for initial 
data in region B. The steady solution was asymmetric over 
_____ 

We have been informed by D. Houghton that S. C. Mehrotra 
has communicated the fact that when the depth of the upper 
layer is much larger than the depth of the lower layer, then 
the size of region A is much smaller. 
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the mountain top. This asymmetry was possible since ■ 0, 
at the top o£ the mountain where H x - 0. A schematic drawing 
of the resulting flow is taken from a figure of 13] 
and reproduced as Figure 3 here. The characteristic curves 
dx/dt * are said to be external for i « l r 4 and internal 
for i * 2,3. A confluence of characteristics with i » 1 
(or with i ■ 4) is said to produce an external shock; 
a region having straight characteristics for i » 2 (or for 
i *3) is said to be an internal rarefaction wave, while if 
i ■ 1 or 4 the wave is called external, on the lee side 
of the mountain , p 3 also approaches zero at a point on the 
mountain side. This caused a confluence of the internal 
characteristics which accompany the stationary shock, that 
is thus called internal and which is located on the mountain 
side. 

Note that in a steady state solution the conservation 
of mass in the lower layer, described by the second component 
of equation (1) , implies that <f>u 5 m - constant. Hence, where 
the depth <f> decreases, the velocity u must increase. Observe 
that this situation occurs on the lee side of the mountain, 
and can explain how the high velocities develop there. 

But we remark that the numerical method of [3] 
produced a spurious solution for data in region B'; namely, 
the initially impulsive motion for B' did not seem to tend to 
a steady state over the mountain. In fact the curve 
separating regions B and B' was determined by numerical 
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experiments. That is, for region B ', numerically spurious 
waves were continuously generated on the upstream side of the 
mountain and these waves traveled both up and downstream, 
thus preventing the development of a steady motion. See typical 
Figures 4a and 4b for computer plots of the heights of the 
two air layers at two different times as found by the method 
of [3]., . On the other hand, when using any of the three 
different numerical methods described in the next section, 
the spurious waves do not appear for data in region B'. 

See Figures 5, 6, and 7 that show the solutions obtained by 
these methods, for the same data in region B'. 

We find that for initial data in region B', the character- 
istic ^2 = ' 0 at x = -a, i.e. at the upstream edge of the 
mountain. (Since the flow is not compressing here, the pseudo- 
viscosity term was not switched on in [3] . But, as is 
well known the Lax-Wendrof f scheme is not dissipative where 
a characteristic is stationary, so that the large force term 

H , caused disturbances to begin at x = - a and move into 

0 

the flow. This error growth could be minimized by using a 
much smaller spatial interval, and could be lessened somewhat 
by introducing a smooth transition region between the flat 
earth and the parabolic mountain edge.) We further observe 
that the flow is symmetric over the top of the mountain, until 
a weak stationary shock forms on the lee mountainside, with 
the accompanying confluence of the characteristics dx/dt = 

A stronger downstream shock is also formed with the confluence 
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of the u 2 characteristics. This stronger shock is stationary 
on the mountain side for some range of initial parameters 
in B*, but in the remainder of the B' region this strong shock 
keeps moving slowly downstream. 

In summary, we find that for any data in regions B and B', 
the wind speed near the ground on the lee side is much higher 
than on the upstream side, in the steady state solution that ' 
rapidly develops. 
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3. Numerical Methods 
(1) Pseudo-Viscosity 

The pseudo-viscosity method [6] devised by 
von Neumann and Richtmyer, was applied to the Lax-Wendroff (L-W) 
scheme in 13 J . Here we use W to represent the discrete 

approximation to the solution w of equation (1) . Symbolically 
we set 

(8) W(t + At) « W(t) + At W fc (t) + 1 -jft w tt (t) ' 

where the time derivatives are to be replaced by spatial 
derivatives obtained from equation (1) , and then the first 
and second order spatial derivatives are to be replaced 
by suitably centered difference approximations. Therefore 

(9) W(t + At) = W - At(G x + K) + L^ 2 - {[A(G X + K)J X + K fc } , 

where A 5 3G/3w is the Jacobian matrix and where dll terms 
on the right-hand side of (9) are to be evaluated at W(t) , with 



m x (r^’+H) x + r4>m’ xx 
0 

i& x ($ + H) x + ^' m xx 


The difference expression for (9) is 


(10) W? +1 - w!j - £ IA 0 g" + 2Ax kJ] 
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where for any function <Kx,t) , we define for integer j, 
♦j 5 <MXj,t n ) , X » At/Ax, x^ * j Ax, t R * n At, with Ax 
and At representing the space and time increments used 
in the difference scheme •, the operators A Q and A + are 
defined by 


A 0 f 


£ j+i ■ f j-i > 


4 + £ j - £ j+i - £ j > 


while 


Ax K j+!a = g 


( ♦ j 4 *V r *j + Bj) 

0 
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with the convention for f = <J>, <p', or A, 


£ J* 5 5 ;£ j + l + V ' 


and define Kj_^ and fj_^ by replacing j by j-1; furthermore, 
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2 Ax K. = | 


( *j+l + ♦j-JL* V r *j + H j> 
0 

<$j+l+ *j-l* A 0^j + H j* 
0 


and finally. 


(Ax) 2 (K fc ) j 5 | 


(Agra^) [A 0 (r^!+H.j) ] + 2r ^ju^j-l* ^V^-l* 

0 

<A 0 m'nAo<*j+ Hj) J + 2(*j| +x + ^. x ) A + (A + m.. 1 

0 


The Lax-Wendroff scheme (9) is modified by the introduc- 
tion of a pseudo-viscosity expression B in the At term as 
follows; 


(11) W(t+At) « W - At KG + B) x + K] + 


(At) 2 

~~1 


UA(G X +K) J x + K t > 


where 



Here (Ax) B x is replaced by (B^- ^) , with B..^ defined 

in a manner analogous to the way Kj+^ are defined. 

In [3] , a M at any mesh point j where either u < 0 

X " 1 ~ 

or < 0; while, if neither inequality holds, then a = 0. 


This is the compression switch, which is tested by checking 
centered velocity differences for each layer. But, as explained 
in section 2, we found that it is desirable to keep a / 0 
everywhere. For the parabolic mountain the optimum range of a 
is (1, 2.5). Note that the introduction of the dimensionless 
parameters in (7) for velocity and length implies that a dimen- 
sionless time unit is also introduced consistently by introduc- 
ing /<f> _a/9 as the unit of time. In. this way, it is easy 
to verify that the dimensionless gravitational constant g * 1. 
Our calculations are performed with a = 1 and with periodic 
boundary conditions imposed at x = + L, where L =* 50 snd Ax = 
The impulsive motion could be computed as long as the 
periodicity condition does not produce interference with 
the flow near the mountain. Occasionally, L = 100 was used 
in order to prevent such interference effects. The 
Courant-Friedrichs-Lewy (CFL) stability condition for (L-W) is 
Mmax |y. |) <1. In the pseudo-viscosity scheme, it is 

j 3 

sometimes necessary to restrict the time step further, in 
order to maintain the stability of the scheme. 

(ii) Filter 

The filter was developed in [1] to facilitate 

the calculation of one and two dimensional aerodynamical flows 
with shocks. The filter may be inserted easily into an exist- 
ing program. For example, if the basic Lax-Wendroff scheme 
(9) , (10) is represented in the operator form 
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( 12 ) 


,n 


wj +1 « T Wjj , 


then the filter scheme has the form 


( 13 ) 


V^ +1 ■ T wj , 


n+1 


(14) Wj +1 - v? +1 + £ ( 0 j + ^(v j+1 - Vj +1 ) - e "_j s (Vj +1 - v?+J)j 


That is. 


«5 +1 ■ v j +l + ? t e " +!s A + ^ +1 - 8 !h 5 \ ^-i 1 • 


Note that if 0 = 1, then the filter is called a Shuman filter, 
and it smooths short wave length oscillations in the solution 
everywhere. Here, we employ 9 as a switch; namely, 


(15) 
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j+*5 
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I'w z j 1 i 
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where z is any dependent function that is a good "sensor" 
of large gradients, for example, z s u the velocity in the 
lower layer; m is one less than the order of accuracy of 
the operator T in (12), that is m = 1; and the constant 3 
lies in the range (0,2). For this range of 3 [1] 

showed that . for constant coefficients, the scheme (13) , (14) 
is stable if (12) is stable. In practice, ws found that for 
the two-layered flow the constant 3 should be chosen in the 
range (0.25, 0.75), with the additional proviso that 3=0 


JL5 


at any times when the flow is so smooth that 


“* l A + u jl “ — J — • 

(iii) Hybrid 

Hybrid schemes were developed and analyzed in [ 2 ] 
to facilitate the calculation of flows with shocks. The idea 
is to use an operator that "interpolates" between two operators: 
a first older accurate operator that provides mono- 
tonic , and slightly diffused shock profiles? and a second 
order accurate operator. Here the interpolation is done with 
a switch analogous to (15) which is sensitive to large 
gradients in the solution. Let the first order Lax-Friedrichs 
scheme be denoted by 

(16) W? +1 = S W? 

where 


S W " = I (»"+! + W "-l> + w t 

" I tW j+l + W j-1> * it(G * + K > j 

- 2 (w j + l + w j-l> - 5 ‘ 4 o G j + 2 K j ! 

“ "j + I (4 + W j " V'j-l > - 3 Ii 0 G j + 2 Ax k!?1 . 


16 


Then the hybrid scheme is suggested by the form 

w” +1 * (0S + (1-0) T)W? . 


That is. 


(17) W? +1 = W” - £ [A q g” + 2 Ax KV] + \ [0 j + ^A + Wj - 

+ r [< 1 -8- +!s )A- +Js (i + G"+4x K^l 




- V s "-! + 4x K j-%> + ( 1 - S j> < K t> 


where 


e ? 5 1 (e j +!5 + "W • 


As shown in [1] and [2] , the filtering and hybrid 
schemes are stable when they use the CFL condition appropriate 
to the basic operators S and T. But, the pseudo-viscosity 
scheme may impose a further restriction on the time step 
(see 16]). In fact, we find that when a strong shock occurs 
in the two-layered model, it is sometimes necessary to reduce 
the time step by a factor of 1/2 , to achieve stability for 
the artificial viscosity scheme. 

In conclusion, we exhibit a calculation of the flow 
over a profile of the Rocky Mountains made with the filtering 
method (14) where the constant 3 = 0.4. The topographical 
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cross-section is the one used in [3] with initial data given 
on p. 45. The constant $ -e0 » 9,000 ft., and represents 
the initial depth of the lower layer which lies above the 
5,000 foot elevation of the plains east of Denver, while 21,000 
feet is the initial depth of the upper layer. The initial 
velocities are u = 18 m/sec , and u' = 28 m/sec , while 
r * 0.946 is the ratio of the densities of the two layers. 
Figure 8 shows the position of the tops of the two layers 
after 3750 time steps, which represents about 5 hours and 
25 minutes of elapsed time. The relatively steady flow now 
exhibits the large increase in velocity in the lower layer 
on the lee side of the mountain range, and a "stationary shock" 
that is near the base of the mountain. The mountain profile 
is the lowest curve drawn. Here the 340 km topography is 
taken for a slice from Denver, at x = 1, west to the Colorado 
border at x = -4. The plains to the east of Denver are 
taken to be level at 5000 ft. above sea level, and the slope 
to the west of the Utah-Colorado border is simulated as being 
constant until the height of 5000 ft. is reached at x = - 6.5, 
and further west the ground remains flat. Zero on the Y axis 
corresponds to 5000 ft. above sea level, while 3.32 on the Y 
axis corresponds to 35,000 ft. ahove sea level. 
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CAPTIONS 

Figure 

X Schematic drawing of a two-layered atmosphere over a bump. 

2 Regions in the (M c ,Fq) plane of initial data, for r * 0.8. 
The regions of physically realistic initial data 

are A, B, and £'. The shallow water model is appropriate 
for data in regions B and B'. 

3 Schematic drawing of a flow that results from initial 
data in region B. The curves plotted show H, <p + H, 
and + $ + H as functions of x. The flow is steady 
near the mountain. 

4a, b Computer plots of solution obtained by the method in [3], 
for initial data, F Q * 0.25, M c * 0.6, r = 0.8, in 
region B ' . Note spurious wavelike motion that spreads 
from upstream edge of mountain (a) at 1500 time steps and 
(b) at 1750 time steps. Here artificial viscosity 
coefficient a = 2, a = 1, Ax * 0.05, and \ < 0.85 (of 
the maximum permissible CFL ratio) . At varies in time. 
Dimensionless time T is printed. 

5 For the initial data and parameter values used in Fig. 4, 
computer plot at 1750 time steps of the solution 
obtained with the pseudo-viscosity method (i.e. artificial 
viscosity term used everywhere) . Flow is steady near 

the mountain. 

6 For the initial data used in Fig. 4 and 8 = 0.5, computer 
plot at 1750 time steps of solution obtained with the 
filtering method. Flow is steady near the mountain. 


21 


For the initial data used in Fig. 4 and 0 » 0.25, 
computer plot at 1750 time steps of solution obtained 
with the hybrid method. Flow is steady near the 
raountain. 

Flow over the Rocky Mountains after 3750 time steps, 
about 5 hours and 25 minutes after the initial time. 
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